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Abstract 

Transfer equation in a geometrically thin accretion disk is reexamined under the plane-parallel approx- 
imation with finite optical depth. Emergent intensity is analytically obtained in the cases with or without 
internal heating. For large or infinite optical depth, the emergent intensity exhibits a usual limb-darkening 
effect, where the intensity linearly changes as a function of the direction cosine. For small optical depth, 
on the other hand, the angle-dependence of the emergent intensity drastically changes. In the case without 
heating but with uniform incident radiation at the disk equator, the emergent intensity becomes isotropic 
for small optical depth. In the case with uniform internal heating, the limb brightening takes place for 
small optical depth. We also emphasize and discuss the limb-darkening effect in an accretion disk for 
several cases. 

Key words: accretion, accretion disks — black holes — galaxies: active — radiative transfer — 
relativity 



1. Introduction 

Accretion disks are now widely believed to be energy 
sources in various active phenomena in the universe: 
in protoplanetary nebulae around young stellar objects, 
in cataclysmic variables and super soft X-ray sources, in 
galactic X-ray binaries and microquasars, and in active 
galaxies and quasars. Accretion-disk models have been 
extensively studies during these three decades (see Kato 
et al. 1998 for a review). 

Radiative transfer in the accretion disk has been inves- 
tigated in relation to the structure of a static disk atmo- 
sphere and the spectral energy distribution from the disk 
surface (e.g., Meyer, Meyer- Hofmeister 1982; Cannizzo, 
Wheeler 1984; Kfiz and Hubeny 1986; Shaviv, Wehrse 
1986; Adam et al. 1988; Hubeny 1990; Ross et al. 1992; 
Artemova et al. 1996; Hubeny, Hubeny 1997, 1998; 
Hubeny et al. 2000, 2001; Davis et al. 2005; Hui et 
al. 2005). In many cases the diffusion approximation 
or Eddington one was employed; it provides a satisfactory 
description at large optical depth, although the emergent 
radiation field originates at optical depth of the order of 
unity. Furthermore, gray and non-gray models of accre- 
tion disks were constructed under numerical treatments 
(Kriz and Hubeny 1986; Shaviv and Wehrse 1986; Adam 
et al. 1988; Ross et al. 1992; Shimura and Takahara 1993; 
Hubeny, Hubeny 1997, 1998; Hubeny et al. 2000, 2001; 
Davis et al. 2005; Hui et al. 2005) and under analytical 
ones (Hubeny 1990; Artemova et al. 1996). 

In these studies, however, the vertical movement and 
the mass loss were not considered. Hence, recently, in re- 
lation to the radiative disk wind, radiative transfer in a 
moving disk atmosphere was also investigated (e.g., Fukue 
2005a, b, 2006a, b). In contrast to the static atmosphere, 
in the moving atmosphere the boundary condition at the 



surface of zero optical depth should be modified (Fukue 
2005a, b). Moreover, the usual Eddington approximation 
violates in the highly relativistic flow (Fukue 2005b; see 
also Turoha, Nobih 1988; Turoha et al. 1995; Dullemond 
1999), and the velocity-dependent variable Eddington fac- 
tor was proposed (Fukue 2006b). 

In the usual studies of radiative transfer in the disk, 
the emergent intensity has not been fully obtained, since 
the attention were usually focused on the disk internal 
structure such as a temperature distribution. In addition, 
the effect of limb darkening has not been well examined 
except for a few cases relating to cataclysmic variables 
(e.g., Diaz et al. 1996; Wade, Hubeny 1998; see also Fukue 
2000; Hui et al. 2005 for high energy cases). 

In this paper, we thus reexamine radiative transfer 
in the accretion disk with finite optical depth under the 
plane-parallel approximation, and analytically obtain an 
emergent intensity for the cases with or without internal 
heating. Besides cataclysmic variables, we also emphasize 
and discuss the limb-darkening effect in an accretion disk 
for various cases. 

In the next section we describe the basic equations. In 
section 3, we show analytical solutions. In section 4, we 
discuss several cases of accretion disk models, and empha- 
size the importance of the limb-darkening effect. The final 
section is devoted to concluding remarks. 

2. Basic Equations 

We here assume the followings: (i) The disk is steady 
and axisymmetric. (ii) It is also geometrically thin and 
plane parallel, (iii) As a closure relation, we use the 
Eddington approximation. (iv) The gray approxima- 
tion, where the opacity does not depend on frequency, 
is adopted, (v) The viscous heating rate is concentrated 
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at the equator or uniform in the vertical direction. 

The radiative transfer equations are given in several 
literatures (Chandrasekhar 1960; Mihalas 1970; Rybicki, 
Lightman 1979; Mihalas, Mihalas 1984; Shu 1991; Kato et 
al. 1998). For the plane-parallel geometry in the vertical 
direction (z)^ the frequency-integrated transfer equation, 
the zeroth moment equation, and the first moment equa- 
tion become, respectively. 
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where 6 is the polar angle, / the frequency-integrated spe- 
cific intensity, E the radiation energy density, F the verti- 
cal component of the radiative fiux, P the 2;z-component 
of the radiation stress tensor, p the gas density, and c the 
speed of light. The mass emissivity j and opacity /^abs 
and A^sca are assumed to be independent of the frequency 
(gray approximation). 

For matter, the vertical momentum balance and energy 
equation are, respectively, 
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where tjj is the gravitational potential, p the gas pressure, 
and the viscous-heating rate. In this paper, we do not 
solve the hydrostatic equilibrium (4). Generally speaking, 
when the contribution of the radiative fiux is small, com- 
pared with the pressure gradient term, the gas pressure 
dominates in the atmoshere, and the density distribution 
will not be constant. When the radiative fiux is strong, 
on the other hand, the radiation pressure dominates, and 
the density may be approximately constant throughout 
much of the disk. Anyway, we suppose that the density 
distribution would be a justed so as to hold the hydro- 
static equilibrium (4) through the main part of the disk 
atmosphere, under the radiative fiux obtained later. 

Using this energy equation (5) and introducing the op- 
tical depth, defined by 

dr = -p (A^abs + i^sca) dz, (6) 
we rewrite the radiative transfer equations: 
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where p = cosO. Final equation is the usual Eddington 
approximation. 



As for the boundary condition at the disk surface of 
r = 0, we impose a usual condition: 



3cP^ = cE^ = 2F^ at r = 0, 



(11) 



where the subscript s denotes the values at the disk sur- 
face. 

For the internal heating, we consider two extreme cases: 
(i) No heating (g^^ = 0), where the viscous heating is 
concentrated at the disk equator and there is no heat- 
ing source in the atmosphere, (ii) Uniform heating in the 
sense that q'^:^^/{Hi^^^ + Hi^'^^)p =constant. The latter case 
means that the kinematic viscosity u is constant in the 
vertical direction, since Q,^-^^/ p = v(rdQ./ dr^ , as long as 
the opacities are constant. 

Finally, the disk total optical depth becomes 



(12) 



ro = - p(/^abs + l^sca)dz, 
JH 

where H is the disk half- thickness. 
3. Analytical Solutions 



Except for the emergent intensity /, several analyti- 
cal expressions for moments as well as temperature distri- 
butions were obtained by several researchers (e.g., Laor, 
Netzer 1989; Hubeny et al. 2005; Artemova et al. 1996). 
For the completeness, we recalculate them as well as the 
intensity /. 

3.1. No Heating Case 

We first consider the case without heating in the disk 
atmosphere: = 0, but with uniform incident intensity 
/o from the disk equator, where the viscous heating is 
assumed to be concentrated. 

In this case, the analytical solutions of moment equa- 
tions are easily given as 



F = K = tt/o 



3cP = c£; = 3Fs I - + r 



(13) 
(14) 



This is a familar solution under the Milne-Eddington ap- 
proximation for a plane-parallel geometry. It should be 
noted that the vertical radiative fiux F is conserved, and 
equals to tt/q at the disk equator. 

Since we obtain the radiation energy density E in the 
explicit form, we can now integrate the radiative transfer 
equation (7). After several partial integrations, we obtain 
both an outward intensity /(r, p) {p > 0) and an inward 
intensity /(r, —p) as 
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where I{ro,p) is the boundary value at the midplane of 
the disk. 
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Fig. 1. Normalized emergent intensity as a function of /i for 
the case without heating. The numbers attached on each 
curve are values of tq at the disk midplane. The dashed 
straight line is for the usual plane-parallel case with infinite 
optical depth. 

In the geometrically thin disk with finite optical depth 
To and uniform incident intensity Iq from the disk equator, 
the boundary value /(ro,/i) of the outward intensity / 
consists of two parts: 

/(ro,/i)=/o + /(ro,-/i), (17) 

where Iq is the uniform incident intensity and I{ro^—fi) is 
the inward intensity from the backside of the disk beyond 
the midplane. Determining I{tq^—ii) from equation (16), 
we finally obtain the outward intensity as 
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where we have used = tt/q . 

For sufficiently large optical depth tq, this equation (18) 
reduces to the usual Milne-Eddington solution: 

-m). (19) 



47r 



Finally, the emergent intensity /(0,/i) emitted from the 
disk surface for the finite optical depth becomes 



-ro//x 
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In figure 1, the emergent intensity l{0^jii) normalized by 
the isotropic value / (= Fs/n) is shown for several values 



of To as a function of /i. 

As is easily seen in figure 1, for large optical depth 
(ro > 10) the angle-dependence of the emergent intensity 
is very close to the case for a usual plane-parallel case with 
infinite optical depth. Therefore, the usual limb-darkening 
effect is seen. Namely, in the case of a semi-infinite disk 
with large optical depth, the energy density increases lin- 
early with the optical depth in the atmosphere, and the 
temperature increases accordingly. As a result, an ob- 
server at a pole-on position of /i = 1 will see deeper in the 
disk, where the temperature (and therefore, the source 
function) is larger than that observed by an observer at 
an edge-on position of /i = 0. Thus, the observed intensity 
will be higher at /i = 1. This is just a usual limb-darkening. 

For small optical depth, however, the angle-dependence 
is drastically changed. When the optical depth is a few, 
the vertical intensity (/i ~ 1) decreases due to the finiteness 
of the optical depth. That is, we cannot see the 'deeper' 
position in the atmosphere, compared with the case of a 
semi-infinite disk. Furthermore, when the optical depth 
is less than unity, the intensity in the direction of small 
/i increases, and the emergent intensity becomes isotropic 
with a uniform value Iq at the disk equator; the limb- 
darkening effect disappears. Indeed, the limiting case of 
To ~ 0, /(0,/i) ~ Fs/tt. That is, in this case of very small 
optical depth, the source function is dominated by the 
isotropic source at the midplane. 

3.2. Uniform Heating Case 

Now, we consider the case with uniform heating: 
<?vis/(^abs + ^sca)p =constant. 

Integrating the equation (8) under the following bound- 
ary conditions: 



F = at r = tq^ 
F = Fs at r = 0, 
we obtain 



F = FJ1 



To 



(21) 



(22) 



The radiative ffux F linearly increases from to the sur- 
face value Fg. 

Substituting equation (22) into equation (9), and inte- 
grating the resultant equation under the boundary condi- 
tion (11), we obtain 

3cP = cE = 3F.(^l+r-Q. (23) 

This expression (23) for finite optical depth is seen in, 
e.g., Laor and Netzer (1989). A similar but more general 
expression was obtained by Hubeny (1990). In any case, 
this expression reduces to the Milne-Eddington solution 
for sufficiently large optical depth. In the case of finite 
optical depth, the radiation energy density and pressure 
decrease from the midplane to the surface in the quadratic 
form. It should be noted that at the midplane of the disk 
of r = To, 
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(24) 
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As already mentioned by Hubeny (1990), the energy den- 
sity at the disk midplane is the half of the corresponding 
stellar atmospheric one. This is explained by the fact that 
the radiation from the disk midplane may escape equally 
to both sides of the disk. 

Since we obtain the radiation energy density E in the 
explicit form (23), we can now integrate the radiative 
transfer equation (7). After several partial integrations, 
we obtain both an outward intensity /(r, fi) {fi > 0) and 
an inward intensity /(r, — /i) as 

^2 



/(r,/i) = 



3^ 
47r 



T 

T 



+/(ro,/i)e 



(r- 



(25) 



T 

~2 



- - 

TO V3 

3to 

-To)/Ai 

1 

To 

3ro To 

where I{ro^fi) is the boundary value at the midplane of 
the disk. 

In the case with uniform heating and without the inci- 
dent intensity, the boundary value /(ro,/i) of the outward 
intensity / is 

/(ro,/i)=/(ro,-/i), (27) 
and we finally obtain the outward intensity as 
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For sufficiently large optical depth ro, this equation (28) 
also reduces to the usual Milne-Eddington solution (19). 

Finally, the emergent intensity /(0,/i) emitted from the 
disk surface for the finite optical depth becomes 
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In figure 2, the emergent intensity /(0,/i) normalized by 
the isotropic value / (= Fq/tt) is shown for several values 
of To as a function of /i. 

As is easily seen in figure 2, for large optical depth 
(ro > 10) the angle-dependence of the emergent intensity 
is very close to the case with a usual plane-parallel case 
with infinite optical depth. Therefore, the usual limb- 
darkening effect is seen, as already stated at the end of 
section 3.1. 

For small optical depth, however, the angle-dependence 




Fig. 2. Normalized emergent intensity as a function of /i 
for the case with uniform heating. The numbers attached on 
each curve are values of ro at the disk midplane. The dashed 
straight line is for the usual plane-parallel case with infinite 
optical depth. 

is drastically changed similar to the case without heat- 
ing. In the vertical direction of /i ~ 1, the emergent in- 
tensity decreases as the optical depth decreases. This is 
due to the finiteness of the optical depth. That is, we 
cannot see the 'deeper' position in the atmosphere, com- 
pared with the case of a semi-infinite disk. In the inclined 
direction of small /i, on the other hand, the emergent in- 
tensity becomes larger than that in the case of the infinite 
optical depth. Moreover, when the optical depth is less 
than unity, the emergent intensity for small /i is greater 
than unity: the limb brightening takes place. Indeed, in 
the limiting case of ro ~ 0, /(0,/i) ~ (Fs/7r)/(2/i). This is 
because that the path length is longer for such a case of 
small /i. That is, in this case for low optical depth, the 
source function is very uniform. This, coupled with the 
absence of an isotropic source at the midplane, is why the 
geometric effect (longer path length) is dominant and one 
finds limb 'brightening'. 

3. 3. Validity of the Eddington Approximation 

In this subsection, we briefly discuss the validity of the 
closure relation in the present treatment. In this paper, we 
have adopted the usual Eddington approximation, where 
the ratio of the radiation pressure to the energy density 
is fixed as 1/3, to close the moment equations. As is well 
known, this approximation is correct in the limit of an 
isotropic radiation field. Hence, in the problem of limb- 
darkening, where the deviation from the isotropy is es- 
sential, this approximation is only approximately correct, 
although it is used in the usual Milne-Eddington approx- 
imation. 

For example, let us suppose the case of a semi-infinite 
disk with an infinite optical depth ro. In this case, we 
easily calculate the energy density as well as the radia- 
tion pressure from the derived intensity (19) of the Milne- 
Eddington solution. At the deeper position in the atmo- 
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sphere, where the integration is done in ah directions, the 
re-calculated variables satisfy the condition oi P/ E = !/?>. 
However, at the surface of the disk, where the integration 
should be done in a semi-sphere, this is not true. At the 
surface of r = 0, we integrate the emergent intensity to 
yield: 
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or P/£^ = (1/3) (17/14). Hence, in the case under the usual 
limb-darkening effect, the ratio of the radiation pressure 
to the energy density is slightly larger than 1/3 at the 
surface. This is just the peaking effect originated from 
the anisotropic radiation field near to the disk surface. 

On the other hand, in the limb-brightening case of small 
optical depth with uniform heating, the situation is re- 
versed. In this case, the ratio of the radiation pressure 
to the energy density becomes smaller than 1/3, and the 
Eddington approximation holds approximately. This is 
also originated from the anisotropy, and may be called 
an anti-peaking effect. As is seen in figure 2, the limb- 
brightening becomes stronger and stronger for small op- 
tical depth. Hence, for such a case of very small optical 
depth, the Eddington approximation would not be good, 
although the qualitative properties would not be changed. 

In order to obtain the intensity distribution more pre- 
cisely, we, for example, introduce a variable Eddington 
factor, that is beyond the scope of the present paper. 

4. Discussion 

As was derived in the previous section, the emergent 
intensity / of the accretion disk depends on the disk total 
optical depth tq as well as the direction cosine fi. The 
limb-darkening effect is considerably modified for small tq, 
compared with the usual case for infinite optical depth. 
Even for the case with sufficiently large optical depth, 
limb darkening in the luminous accretion disk must be 
important, and should be examined more carefully. 

In this section, we discuss several cases in turn, and 
call the attention to the importance of the limb-darkening 
effect. 

4.1. Standard Disk 

The optical depth at the midplane in the inner region 
of a geometrically thin standard disk (Shakura, Sunyaev 
1973; see also Kato et al. 1998) is expressed as 



To ^A^i; : 
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where k is the electron scattering optical depth, a the 
viscous parameter, m the mass accretion rate normalized 
by the critical rate M (=Le/c^), i^E being the Eddington 
luminosity of the central object, f the radius normalized 
by the Schwarzschild radius rg (= 2GM/c^). 



This optical depth becomes small for large m and/or 
small f. For a slightly large accretion rate, inside some 
critical radius 

Tcr = 2m, (33) 

the disk shifts to a supercritical regime, whereas the disk 
is a standard regime outside r^r (Fukue 2004). At this 
critical radius, the optical depth becomes 

Tcr - 57a-^m^/^ (34) 

Hence, in a usual situation the optical depth of the inner 
region of the standard disk is greater than several tens. 

In such a situation, however, due to a usual limb- 
darkening effect for semi-infinite medium, the emergent 
radiation toward the pole-on direction is enchanced by 
20 percent, while the emergent radiation seen from the 
edge-on direction diminishes by 50 percent. Thus, in cal- 
culating the flux and spectrum of the standard disk, we 
carefully consider the limb-darkening effect. 

In the region inside the inner edge at 3rg, the disk gas 
freely falls toward the central black hole, and the surface 
density (i.e., the disk optical depth) quickly drops. Hence, 
the emergent spectrum from the innermost region inside 
the inner edge would be greatly modified from the opti- 
cally thick case and the optically thin one. 

4-2. Supercritical Disk 

In the supercritical accretion disk, where the mass ac- 
cretion rate exceeds the critical rate, the expression for the 
disk optical depth is changed. For example, in the self- 
similar model without mass loss (Fukue 2000), the optical 
depth at the midplane of the disk is 



n M 1 



m 1 



'^o = -; , = —r -5 (35) 

where ci is a coefficient of the order of unity. At the 
critical radius, the disk optical depth is 



1 



Tcr ~ ■ 



2cia 



-m 



1/2 



(36) 



In the critical accretion disk (Fukue 2004), where the 
mass accretion rate exceeds the critical one, but the excess 
mass is expelled by the wind mass loss, the optical depth 
at the midplane of the disk is 



To 



16a/6.i 



pi/2 = 39.2a-ifi/^ 



(37) 



At the critical radius, the disk optical depth is 

rcr-55a-^m^/^ (38) 

Hence, in a usual situation the optical depth of the su- 
percritical/critical disk is also greater than several tens. 

In such a situation, however, the usual limb-darkening 
effect for semi-infinite medium is also important (see 
Fukue 2000). In addition to the limb-darkening effect, 
the geometrical effects, such as a projection effect and a 
self-occultation, should be considered in calculating the 
flux and spectrum of the supercritical disk (e.g., Fukue 
2000; Watarai et al. 2005; Kawata et al. 2006). 
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4.3. Disk Corona and ADAF 

If a luminous disk is sandwiched by a disk corona, the 
situation is similar to the case without heating, but with 
the incident intensity from the midplane. Hence, when 
the optical depth of the corona is sufficiently smaller than 
unity, the emergent intensity is just that of the disk, ex- 
cept for a very small direction cosine. However, when the 
optical depth of the corona is a few, limb darkening takes 
place, and the emergent intensity toward the edge-on di- 
rection remarkably reduces. 

On the other hand, if the accretion rate is quite small, 
and the inner region of the disk becomes an optically- 
thin advection dominated state (ADAF), the situation is 
similar to the case with internal heating, although, rig- 
orously speaking, the plane-parallel approximation may 
be invalid. When the optical depth of an ADAF region is 
sufficiently smaller than unity, the limb brightening would 
take place. 

In both cases with hot gas, however, the effect of the 
Compton scattering may be important. Hence, the angle 
dependence of the intensity would be modified by the an- 
gle dependence of the Compton scattering, and the trans- 
fer problem should be treated more carefully. 

In addition, in the latter case of ADAF, the accretion 
flow is supposed to be conical or spherical. Hence, when 
the opening angle is small, the limb brightening would 
qualitatively take place. When the opening angle is large, 
on the other hand, the extension of the emitting region is 
large, and the angle dependence of the emergent intensity 
would become much more complicated than the present 
simple case. 

4^4' Relativistic Disk 

In the case of the relativistic standard disk (e.g., 
Novikov, Thorne 1973; Page, Thorne 1974), the situa- 
tion is similar to the non-relativistic case. However, the 
direction cosine of the local emergent intensity in the disk 
to the observer is changed by two additional reasons: (i) 
the light trajectory is bent by the space- time curvature, 
and (ii) the emission from the gas rotating around a black 
hole suffers from the special relativistic aberration. Limb- 
darkening effect in the comoving frame on the spectral 
energy distribution (SED) from the relativistic accretion 
disk was considered by, e.g., Fu and Taam (1990) and 
Gierlihski et al. (2001). 

Recently, due to submilliarcsecond astrometry, imaging 
of a black-hole silhouette is expected to become possible in 
the near future at infrared and submillimetre wavelengths. 
The "photographs" of relativistic accretion disks around 
a black hole have been obtained by many researchers 
(Luminet 1979; Fukue, Yokoyama 1988; Karas et al. 1992; 
Jaroszyhski et al. 1992; Fanton et al. 1997; Fukue 2003; 
Takahashi 2004, 2005). In these studies, however, the 
usual limb darkening was not considered. 

Examples of silhouettes of a dressed black hole are 
shown in figures 3 and 4. In figure 3 an edge-on view 
with inclination angle of 80° is shown, while a pole-on 
view is expressed in figure 4. In both figures, the left pan- 




Fig. 3. Edge-on view of a dressed black hole without (left 
panel) and with (right panel) limb darkening. The inclination 
angle is 80°. 




Fig. 4. Pole-on view of a dressed black hole without (left 
panel) and with (right panel) limb darkening. The inclination 
angle is 0°. 

els are for the case without limb darkening, whereas the 
right panels are for the case with limb darkening. 

In the edge-on view (figure 3), the image of a limb- 
darkening disk darkens as expected. This is due mainly to 
the usual limb-darkening effect for small direction cosine. 
Surprisingly, on the other hand, in the pole-on view (figure 
4), the image of a limb-darkeing disk also darkens! This is 
because that the local direction cosine becomes small due 
to the light aberration associated with the disk rotation. 

Thus, in calculating spectra and observed fluxed of rel- 
ativistic disks as well as in taking black hole silhouettes, 
we should carefully consider limb darkening. For example, 
besides Fu and Taam (1990) and Gierlihski et al. (2001), 
Hui et al. (2005) discussed limb darkening in some de- 
tails in their non-LTE calculation of accretion disk spectra 
around intermediate- mass black holes. 

5. Concluding Remarks 

In this paper we analytically solve the radiative trans- 
fer problem of a geometrically thin accretion disk with 
finite optical depth, and obtain analytical expressions for 
emergent intensity from the surface of the disk. For small 
optical depth, the angle-dependence of the emergent in- 
tensity drastically changes from the case with large optical 
depth. In the case without heating but with uniform inci- 
dent radiation at the disk equator, the emergent intensity 
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becomes isotropic for small optical depth. In the case with 
uniform internal heating, the limb brightening takes place 
for small optical depth. We also emphasize the impor- 
tance of limb darkening in the accretion disk study, and 
discuss several cases, including reiativistic silhouettes. 

In the calculation of spectra from, e.g., cataclysmic vari- 
ables, the effect of limb darkening was considered (e.g., 
Diaz et al. 1996; Wade, Hubeny 1998). In the cases of 
high energy and reiativistic regimes, we should also con- 
sider limb darkening (cf. Hui et al. 2005). 

In addition, if there exist intense radiation sources, such 
as neutron stars or radiating jets, irradiation takes place 
and the outer boundary condition changes (e.g., Hubeny 
1990). In such cases under strong incident radiation, limb 
brightening may occur (cf. Stibbs 1971), even for large 
optical depth. 

Finally, if there exists mass loss from the disk surface, 
the radiative transfer problem becomes much more com- 
plicated. In order to examine spectra, fluxes, eclipsing 
light curves, we must calculate the emergent intensity an- 
alytically or numerically. 

This work has been supported in part by a Grant-in- Aid 
for Scientific Research (18540240 J.F.) of the Ministry of 
Education, Culture, Sports, Science and Technology. 
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